Exploring some ode models
# load key packages
library(Rmetapop)
source("../../R/ancillary.R")
source("../../R/survival_straying_recruitment.R")
source("../../R/stochasticity.R")
library(dplyr)
library(R2admb)
library(parallel)
library(R.utils)
library(MASS)
library(VGAM)
library(RColorBrewer)Helper functions
There are a variety of helper functions in this code. Check out some of the helper files:
Spatial correlation
The basic spatial covariance kernel (or correlation where spat_sd = 1) is a squared exponential or periodic.
?spat_cor_mat
example(spat_cor_mat)##
## spt_c_> #a correlation matrix with 10 sites and sd of 4
## spt_c_> round(spat_cor_mat(n_loc=10,spat_sd=4,sumto1=FALSE),3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898
## [2,] 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871
## [3,] 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838
## [4,] 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812
## [5,] 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802
## [6,] 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812
## [7,] 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838
## [8,] 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871
## [9,] 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898
## [10,] 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000
##
## spt_c_> #a stray th 10 sites and sd of 4
## spt_c_> round(spat_cor_mat(n_loc=10,spat_sd=4,sumto1=TRUE),3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104
## [2,] 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101
## [3,] 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097
## [4,] 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094
## [5,] 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093
## [6,] 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094
## [7,] 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097
## [8,] 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101
## [9,] 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104
## [10,] 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116
##
## spt_c_> # illustrate different parameter values
## spt_c_> library(RColorBrewer)
##
## spt_c_> a <- data.frame(n_loc= rep(10,3),spat_scale=c(1,4,10),spat_sd=c(1,4,10))
##
## spt_c_> b <- array(cbind(with(a,mapply(spat_cor_mat,spat_scale=spat_scale,n_loc=n_loc)),
## spt_c_+ with(a,mapply(spat_cor_mat,spat_sd=spat_sd,n_loc=n_loc))),dim= c(10,10,6))
##
## spt_c_> df1 <- melt(b,varnames = c("site1","site2","parameter"),value.name="correlation")
##
## spt_c_> df1$site1 <- factor(df1$site1)
##
## spt_c_> df1$parameter <- factor(df1$parameter, labels=
## spt_c_+ c(paste(rep(c("Cauchy, scale=","Gaussian, sd="),each= 3),rep(c(1,4,10),2))))
##
## spt_c_> ggplot(aes(site2,correlation),data =df1)+
## spt_c_+ geom_line(aes(colour=site1))+
## spt_c_+ labs(colour="site")+
## spt_c_+ scale_x_continuous(breaks= c(1:10))+
## spt_c_+ xlab("site")+
## spt_c_+ facet_wrap(~parameter)+theme_bw()
## spatiotemporal correlation
The basic genesis of a spatiotemporal time series (matrix n time points long for s locations) uses cholesky decomposition (L) of the separable temporal AR(1) and a spatial periodic (circular) covariance functions and a random white noise matrix \(\epsilon\):
\(L_t [L_s \epsilon']'\)
?spat_temp_ts
example(spat_temp_ts)##
## spt_t_> ### generate the time series
## spt_t_> errors <- spat_temp_ts(n_iter = 100, n_loc = 5,
## spt_t_+ site_sd = 0.5, spat_sd = 2,
## spt_t_+ phi = 0.5)
##
## spt_t_> # examine correlation matrix of time series
## spt_t_> errors$cor_mat
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.8260514 0.5969475 0.6017253 0.7671826
## [2,] 0.8260514 1.0000000 0.8163391 0.6186857 0.5632422
## [3,] 0.5969475 0.8163391 1.0000000 0.7765027 0.5311337
## [4,] 0.6017253 0.6186857 0.7765027 1.0000000 0.7644076
## [5,] 0.7671826 0.5632422 0.5311337 0.7644076 1.0000000
##
## spt_t_> # examine pacf of each time series
## spt_t_> errors$pacf
## [1] 0.4938427 0.5033885 0.5299877 0.5191227 0.4103612
##
## spt_t_> # examine sd of each time series
## spt_t_> apply(errors$ts,2,sd)
## [1] 0.5630387 0.5339049 0.5269409 0.5730403 0.4730050
##
## spt_t_> # plot the time series
## spt_t_> matplot(errors$ts,type= "l")
stray matrix
The migration matrix is a random Dirichlet draw from a mean probability and random scale
?ran_stray_prob
example(ran_stray_prob)##
## rn_st_> stray_mat <- spat_cor_mat(n_loc=10,spat_sd=2, sumto1=TRUE)
##
## rn_st_> stray_mat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.13732540 0.11902069 0.10503539 0.08999813 0.07942307 0.07572006
## [2,] 0.11902069 0.13732540 0.11902069 0.10503539 0.08999813 0.07942307
## [3,] 0.10503539 0.11902069 0.13732540 0.11902069 0.10503539 0.08999813
## [4,] 0.08999813 0.10503539 0.11902069 0.13732540 0.11902069 0.10503539
## [5,] 0.07942307 0.08999813 0.10503539 0.11902069 0.13732540 0.11902069
## [6,] 0.07572006 0.07942307 0.08999813 0.10503539 0.11902069 0.13732540
## [7,] 0.07942307 0.07572006 0.07942307 0.08999813 0.10503539 0.11902069
## [8,] 0.08999813 0.07942307 0.07572006 0.07942307 0.08999813 0.10503539
## [9,] 0.10503539 0.08999813 0.07942307 0.07572006 0.07942307 0.08999813
## [10,] 0.11902069 0.10503539 0.08999813 0.07942307 0.07572006 0.07942307
## [,7] [,8] [,9] [,10]
## [1,] 0.07942307 0.08999813 0.10503539 0.11902069
## [2,] 0.07572006 0.07942307 0.08999813 0.10503539
## [3,] 0.07942307 0.07572006 0.07942307 0.08999813
## [4,] 0.08999813 0.07942307 0.07572006 0.07942307
## [5,] 0.10503539 0.08999813 0.07942307 0.07572006
## [6,] 0.11902069 0.10503539 0.08999813 0.07942307
## [7,] 0.13732540 0.11902069 0.10503539 0.08999813
## [8,] 0.11902069 0.13732540 0.11902069 0.10503539
## [9,] 0.10503539 0.11902069 0.13732540 0.11902069
## [10,] 0.08999813 0.10503539 0.11902069 0.13732540
##
## rn_st_> rand_strays1 <- ran_stray_prob(stray_mat,n_iter=3,scale=100)
##
## rn_st_> rand_strays2 <- ran_stray_prob(stray_mat,n_iter=3,scale=1000)
##
## rn_st_> rand_strays3 <- ran_stray_prob(stray_mat,n_iter=3,scale=10^5)
##
## rn_st_> df1 <- melt(rand_strays1,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df2 <- melt(rand_strays2,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df3 <- melt(rand_strays3,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df4 <- rbind(df1,df2,df3)
##
## rn_st_> df4$precision <- factor(rep(c(1,2,3),each=nrow(df4)/3), labels= c("Precision = 10^2", "Precision =10^3", "Precision =10^5"))
##
## rn_st_> df4$site1 <- factor(df4$site1)
##
## rn_st_> ggplot(aes(site2,stray_probability),data =df4)+
## rn_st_+ geom_line(aes(colour=site1))+
## rn_st_+ labs(colour="source site")+ylab("straying probability (random Dirichlet samples")+
## rn_st_+ scale_x_continuous(breaks= c(1:10))+
## rn_st_+ xlab("destination site")+
## rn_st_+ facet_grid(precision~replicate)+theme_bw()
Set parameters for the example
Stocks and substocks
Run a simulation where there are 10 stages (1 egg stage, a 1 year old recruit stage, and ages 2-10). Lets say there is one stock, with six sub-stocks at unique locations.
### basic spatial properties
n_stages= 10 ### number of stages (eggs, recruits, adults)
n_stocks <- 1 ### number of stocks to simulate
n_loc_stocks <- 6 ### substocks per stock
n_loc = n_stocks*n_loc_stocks ### total number of "stocklets"Life history
The maturity of ages 2 and 3 are 25% and 90% followed by 100% thereafter. We can take beverton-holt parameters from the pacific herring stock assessment. We also assume instantaneous natural mortality is 0.58. Assume desnity dependence in recruitment is determined locally, and fish do not follow other fish to spawn (e.g., like MacCall et al. 2018) but this could be done by setting follow to TRUE. Assume no difference in productivity among locations.
maturity = c(0.25,0.9,1,1,1,1,1,1,1) # proportional maturity for ages 2 through plus group
stage_mat <- 2 # stage at first maturity
### taken from Table 13 of the 2015 herring assessment
B0 =26.81/n_loc_stocks # tons in central coast in 2013/14
E0 = B0*1e6*100*1000 # eggs at 6 years old
h_assess <- 0.8 # group steepness
### taken from table E.7 of the 2014 herring assessment
R0 = 287.17e06/n_loc_stocks # maximum local recruitment
beta_0= 1/R0
alpha_0= (E0)/(R0)*(1-h_assess)/(4*h_assess)
dd <- c(1) # density dependence is local
follow <- c(FALSE) # fish do not follow each other (GWOF)
prod_diff = 0.1 # no productivity differential
M_assess <- c(0.58) # natural mortality (fixed)#### stochasticity parmameters
rec_sd <- c(0.8) # log-scale sd in rec-devs
C <- 1e10 ### dirichlet stochasticity in migration (high C is low stochasticity)
obs_sd = c(0.4) ### observation error
phi <- c(0.6) ### recruitment autocorrelation
surv_rho <-1e-6 ### survival stochasticity (low rho is low stochasticity)Simulation parameters
### simulation parameters
n_iter = 50 #length of the time series
start = 21 #when to start assessment
warmup = 20 #
warmup2 = 25 # when to start recording data
n_reps= 100 # number of times to run (not run)
s_timeout= 60 # seconds to timeout and restart for a given assessment run
### assessment parameters (bow many samples to take from the simulated population)
ns_samps=200
nf_samps=200Test functions to esnure they do what you think they’re doing:
### set cauchy scale & recr. corrs to achieve desired retention & syncrony
retention <- function(x)
sapply(x,function(y) mean(diag(spat_cor_mat(n_loc, spat_scale =y, sumto1 = TRUE,offset= 0))))
obj_stray <- 1-c(0.01,0.2,0.4)
obj_juv_stray <- 1-c(0.2,0.4)
ret_fun <- function(x) exp(optim(1,
function(p) (retention(exp(p))-x)^2,
method="BFGS")$par)
stray_s_scale <- sapply(obj_stray,ret_fun)
retention(stray_s_scale)## [1] 0.9900004 0.8000003 0.5999998juv_stray_s_scale<- sapply(obj_juv_stray,ret_fun)
retention(juv_stray_s_scale)## [1] 0.8000003 0.5999998offset=0.1
synchrony <- function(n_loc,spat_sd)
mapply(function(y,x) {(1+(y-1)*(mean(spat_cor_mat(y, spat_sd =x)[lower.tri(spat_cor_mat(y, spat_sd =x))])))/y}, n_loc,spat_sd)
obj_synchrony <-c(0.2)
sd_fun <- function(x) exp(optim(2,
function(p) (synchrony(n_loc,exp(p))-x)^2,
method="BFGS")$par)
rec_corr_sd <- sapply(obj_synchrony,sd_fun)
### pre-check for positive definite covariance matrix using cholesky decompose
lapply(rec_corr_sd,function(x) spat_cor_mat(n_loc,spat_sd= x,offset= offset))## [[1]]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000000 0.0987701771 0.0011659042 0.0001266722 0.0011659042
## [2,] 0.0987701771 1.0000000000 0.0987701771 0.0011659042 0.0001266722
## [3,] 0.0011659042 0.0987701771 1.0000000000 0.0987701771 0.0011659042
## [4,] 0.0001266722 0.0011659042 0.0987701771 1.0000000000 0.0987701771
## [5,] 0.0011659042 0.0001266722 0.0011659042 0.0987701771 1.0000000000
## [6,] 0.0987701771 0.0011659042 0.0001266722 0.0011659042 0.0987701771
## [,6]
## [1,] 0.0987701771
## [2,] 0.0011659042
## [3,] 0.0001266722
## [4,] 0.0011659042
## [5,] 0.0987701771
## [6,] 1.0000000000lapply(rec_corr_sd,function(x) chol(spat_cor_mat(n_loc,spat_sd=x,offset=offset)))## [[1]]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0.09877018 0.001165904 0.0001266722 1.165904e-03 0.0987701771
## [2,] 0 0.99511027 0.099139787 0.0011590603 1.157224e-05 -0.0086318511
## [3,] 0 0.00000000 0.995072833 0.0991436175 1.169158e-03 0.0008715699
## [4,] 0 0.00000000 0.000000000 0.9950724515 9.914263e-02 0.0010823202
## [5,] 0 0.00000000 0.000000000 0.0000000000 9.950719e-01 0.0990348547
## [6,] 0 0.00000000 0.000000000 0.0000000000 0.000000e+00 0.9901313599synchrony(n_loc,rec_corr_sd)## [1] 0.1999998Specify fishing parameters:
### range of fishery parameters for the simulation
Fmort <- c(0,0.1,0.2,0.25,0.375) # fishing harvest rate (proportional)
h_crit= c(0.25) # lower biomass threshold for closure
spat_alloc <- c(1) # spatial allocation options (1 mortality rate is equal across open locations)
stage_forecast= c(2) # first stage included in the forecast
assessment= TRUE # whether or not to use a formal assessmentcombine parameters
Now combine parameter combinations into a dataframe using expand.grid.
### generate the dataframe of parameters and output
param.df <- expand.grid(stray_s_scale= stray_s_scale,
juv_stray_s_scale= juv_stray_s_scale,
rec_corr_sd=rec_corr_sd,
surv_rho= surv_rho,
phi=phi,
Fmort= Fmort,
h_crit= h_crit,
dd= dd,
spat_alloc=spat_alloc,
M_assess=M_assess,
follow= follow,
h= h_assess,
obs_sd= obs_sd,
stage_forecast= stage_forecast,
MPA= c(0,6),
diff= prod_diff)%>%
mutate(a_bh_assess = alpha_0,
b_bh_assess = beta_0,
retention = retention(stray_s_scale),
juv_retention= retention(juv_stray_s_scale),
sync = synchrony(n_loc,rec_corr_sd))
param.df <- subset(param.df,!(follow==TRUE&dd==0)&!(diff!=1&dd==0)&
!(Fmort%in%c(0.025,0.05))&
!(follow==FALSE&juv_stray_s_scale%in%juv_stray_s_scale[c(2,3)])&
!(Fmort==0&h_crit>0.25)&
!(Fmort==0&MPA==6))
param.df[param.df$Fmort==0,"assess"]## NULL# for sd = 0.7, 6 sites
coefs <- c(1.1180355,-1.0114039,0.3912513,1.4970863,-0.9118604,-0.4836503)
param.df$site_sd <- as.numeric(cbind(1,poly(log10(param.df$rec_corr_sd),5,raw= TRUE))%*%t(t(coefs)))
# set ids by row for later use
param.df$id <- 1:nrow(param.df)Run the simulation
Create a wrapper for the simulation function
We’re using an ADMB assessment, but could avoid an assessment using assessment = FALSE.
### give current directory so we can come back after running admb
wd <- getwd()
#### example simulation function ####
example_sim<- function(x){
setwd(wd)
directory= tempfile(tmpdir= getwd())
dir.create(directory)
file.copy(c("Herring.tpl","Herring","Herring.cpp","Herring.obj","Herring.htp"),directory)
setwd(directory)
rep <- fishery_simulate(n_loc=n_loc, ### total number of subpopulations
n_stages=n_stages, ### number of age/stage classes
stage_mat=2, ### stage/age at maturity
maturity = maturity,
a_bh=1/c(rep(1/x$a_bh_assess*(1-(x$diff-1)/n_loc),5),1/x$a_bh_assess*x$diff), ### BH alpha parameter
b_bh=1/c(rep(1/x$b_bh_assess*(1-(x$diff-1)/n_loc),5),1/x$b_bh_assess*x$diff),
phi=x$phi, ### lag 1 recruitment autocorrelation
M=x$M_assess, ### natural mortality (a constant, vector or matrix)
surv_rho= x$surv_rho, ### random survival stochasticity
assessment= TRUE, ### whether to assess the stock
spat_alloc= 2,
collective_dd=x$dd, ### group density dependence
spat_scale=x$stray_s_scale, ### spatial Cauchy straying scale
juv_spat_scale=x$juv_stray_s_scale, ### spatial Cauchy straying scale
spat_sd=x$rec_corr_sd, ### spatial SD of recruitment synchrony
site_sd=x$site_sd, ### degree of stochasticity in alpha
C=C, ### inverse of degree of stochasticity in straying
obs_sd=x$obs_sd, ### observation error on surveys of total biomass
stock_IDs= rep(1:n_stocks,each= n_loc/n_stocks), ### ids representing to which stock the populations belong
harvest_IDs= c(rep(1,n_loc-1),ifelse(x$MPA==0,1,NA)),
n_iter=n_iter, ### total number of iterations
warmup= warmup,
warmup2= warmup2, ### warmup prior to harvest
point.estimate=FALSE, ### assessment point estimate vs full posterior
h_crit=x$h_crit, ### lower harvest cutoff
Fmort=x$Fmort, ### proportional fishing rate
fit.compile=fit.compile, ### compiled assessment model
ret_ts= TRUE,
follow=x$follow,
stage_forecast= x$stage_forecast)
setwd(wd)
unlink(directory,recursive= TRUE)
return(list(TS= rep$TS,metrics = rep$metrics,rep$age_quantiles))
}Run a single simulation
data <- rows.to.list(param.df)
out <- example_sim(data[[20]])## iterations 1:20 of 50 complete (warmup)
## iteration 22 of 50 completed 0.257 sec
## iteration 23 of 50 completed 0.224 sec
## iteration 24 of 50 completed 0.242 sec
## iteration 25 of 50 completed 0.234 sec
## iteration 26 of 50 completed 0.264 sec
## iteration 27 of 50 completed 0.273 sec
## iteration 28 of 50 completed 0.283 sec
## iteration 29 of 50 completed 0.3 sec
## iteration 30 of 50 completed 0.306 sec
## iteration 31 of 50 completed 0.284 sec
## iteration 32 of 50 completed 0.312 sec
## iteration 33 of 50 completed 0.294 sec
## iteration 34 of 50 completed 0.262 sec
## iteration 35 of 50 completed 0.3 sec
## iteration 36 of 50 completed 0.316 sec
## iteration 37 of 50 completed 0.346 sec
## iteration 38 of 50 completed 0.346 sec
## iteration 39 of 50 completed 0.377 sec
## iteration 40 of 50 completed 0.366 sec
## iteration 41 of 50 completed 0.398 sec
## iteration 42 of 50 completed 0.432 sec
## iteration 43 of 50 completed 0.417 sec
## iteration 44 of 50 completed 0.511 sec
## iteration 45 of 50 completed 0.506 sec
## iteration 46 of 50 completed 0.448 sec
## iteration 47 of 50 completed 0.47 sec
## iteration 48 of 50 completed 0.5 sec
## iteration 49 of 50 completed 0.521 sec
## iteration 50 of 50 completed 0.537 secFormat output for plotting
### format example output
index <-reshape2::melt(out$TS$B_stocks)
names(index) <- c("stock","year","type","index_real")
assess <- melt(out$TS$assess)%>%filter(Var3>warmup&Var4=="Estimated Stock Size")
names(assess) <- c("stock","year","iter","type","est")
assess$type <- "retrospective estimate"
forecast <- melt(out$TS$forecast)%>%filter(!is.na(value)&Var4=="Prediction")
forecast$type <- "forecast (pre-harvest)"
sites <-melt(t(out$TS$B[,,1]),value.name= "biomass")%>%
join(melt(out$TS$Loc_Harvest,value.name= "harvest"))%>%
mutate(Site= factor(Var2,labels= c(1,2,3,4,5,"Closed Area")))
site_mean <- sites%>%
group_by(Var1)%>%
dplyr::summarize(mu= sum(biomass),
muH= sum(harvest))%>%
data.frame()
site_mean$B0_est <- out$TS$B0_estGenerate figures
### plot example sims
options <- theme(axis.text.x= element_blank(),
panel.background= element_blank(),
strip.background= element_blank(),
panel.grid.minor= element_blank(),
panel.grid.major.x=element_blank(),
legend.position= c(0.8,0.8),
legend.spacing.y = unit(-0.5, "cm"),
legend.title= element_blank(),
legend.key= element_blank(),
legend.background= element_blank(),
legend.text= element_text(size= 9),
legend.direction= "horizontal",
axis.title= element_text(size= 12),
axis.text.y= element_text(size= 12))
### stock scale plot
stock_plot<- ggplot(aes(year, index_real/out$TS$B0),data=index)+
geom_line(aes(linetype= type,colour= type), size= 1,data= subset(index,type=="actual"))+
geom_line(aes(year,est/out$TS$B0,group= factor(iter),linetype=type,colour= type),
data = subset(assess,type== "retrospective estimate"))+
annotate(x= 10,y= Inf, vjust=1 ,label ="A) stock scale",geom="text")+
geom_segment(aes(y= 0,yend= muH/out$metrics$s_b0
,x= Var1-1,xend= Var1-1),size= 2,col= "grey40",data= site_mean)+
geom_point(aes(Var3,value/out$TS$B0,fill= type),shape= 21,colour= "red",data= forecast)+
geom_point(aes(shape= type),size=1,fill= "grey40",data= subset(index,type== "observed"))+
theme_bw()+
scale_linetype_manual(values= c(1,1),name= "")+
scale_shape_manual(values= c(1),name= "")+
scale_fill_manual(values= c("red","grey40"),name= "")+
scale_colour_manual(values= c("blue",alpha("black",0.4)),name= "")+
geom_line(aes(y= 0.4*B0_est/out$TS$B0,x= Var1-1),linetype= "dotted",col= "red",data= site_mean)+
geom_hline(yintercept= 0.4,col= "red",size= 1,alpha= 0.5)+
ylab("Spawning Biomass")+
#coord_cartesian(ylim= c(0,2),expand= c(0,0),xlim= c(0,n_iter))+
xlab(expression(paste("Time ",symbol('\256'))))+
scale_y_continuous(breaks= c(0,1),labels= c(0,"B0"))+
options+
theme(axis.title.x= element_blank())
### substock scale plot
ages <- t(apply(out$TS$ages,c(1,3),sum))
plot(I(ages[3:50,2]/1e6)~ages[1:48,1],xlim= c(0,max(ages[1:48,1])))
substock_plot <- ggplot(aes(Var1,biomass/out$metrics$p_b0
),data= sites)+
geom_segment(aes(y= 0,yend= muH/out$metrics$s_b0
,x= Var1-1,xend= Var1-1),size= 2,col= "grey40",data= site_mean)+
geom_path(aes(colour= Site,size= Site),lineend= "round")+#facet_grid(Site~.)+
geom_hline(yintercept= 0.25,col= "red",size= 1,alpha= 0.5)+theme_bw()+
ylab("Spawning Biomass")+
annotate(x= 10,y= Inf, vjust=1, label ="B) substock scale",geom="text")+
# coord_cartesian(ylim= c(0,2),expand= c(0,0),xlim= c(0,n_iter))+
xlab(expression(paste("Time ",symbol('\256'))))+
scale_y_continuous(breaks= c(0,1),labels= c(0,"B0"))+
options+
scale_colour_manual(values= c(brewer.pal(7,"PRGn")[c(1:3,6,7)],"black"))+
scale_size_manual(values= c(rep(0.5,5),1))
grid.draw(rbind(ggplotGrob(stock_plot), ggplotGrob(substock_plot), size = "first"))